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SECTION 0. 


SUMMARY 

This is the- written -portion of work performed uncer NASA 
Contract NAS2-5233 - "Report on Pioneer's 6/7 Plasma Probe Instrument 
Simulation Study." In response to the direction of the Technical 
Monitor, vthis report concentrates on part A. 1.(2) of the work state- 
ment "recommendations concerning the design and construction of a second 
generation version of this algorithm ..." The recoirenendati ons contained 
herein are designed to be applicable to the analysis of data from the 
later 3 Collector Electrostatic Analyzer System. As per the contract, 
a minimum of six' days were spent in consultation with Ames Scientists 
during the summer of 1970. 

The results of three independent studies are reported here. They 
are 1.) Optimization of Numerical Quadratures in the Plasma Probe 
Non-Linear Least Squares Parameter Estimation Program, 2.) An 
Alternative Technique for Estimating the Helium Fraction and 3.) 

Linear Estimates of the Plasma Parameters. 
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SECTION 1 . 


INTRODUCTION 

Effort on this contract has been concentrated on the develop- 
ment of techniques for estimating as many as possible of the (steady- 
state) plasma parameters, from data observed by Pioneer 9 type instru- 
ments. The basic technique for doing this for Pioneer 6 and 7 has 
involved a mathematical model of the instrument response, and a math- 
ematical model of the solar wind particle velocity vector distribution, 
as a function of its parameters. A non-linear least squares algorithm 
in conjunction with these models, has been successfully used to esti- 
mate the parameters. Work under this contract has proceeded under the 
assumption that this basic technique is sound for 9 type instruments. 
The 9 instrument response' functions are of somewhat different shape 
than those for Pioneer 6 and 7. 

Loosely speaking, the k th data point 1^ is modeled as 

= M 0 )* where 

V 0 ) “ jfiJ w |<(v)f(v;0)dv » (1.1) 

v is a velocity vector, w^(v) is the total instrument response to a 

unit number density cold beam with velocity vector v at the k^ 1 
choice of azimuth, collector, and energy step combination, and 
f( v-,0) is the particle velocity vector distribution, as a function of 
its m parameters, e = (e^ jQg,.. . ,6 m ). These parameters include the 
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components of the mean vector(s), and pressure tensor(s), number 
density, helium fraction, etc. w^v) is obtained from the calibration 

data and the impulse response of the electronics. 

The non-linear least squares program attempts to iteratively 
choose the parameter set 8 = (9-j, so that the modeled data 

points I|,(8) are close to the corresponding observed data points 1^, 

in an appropriate least squares sense. In order to do this, it is 
necessary to have reasonable starting guesses for the parameters, and 
to evaluate 1^(0) for various parameter sets 0 to some appropriate 

level of accuracy, in an efficient manner. 

The main body of this report consists of three sections. 

Section 2 is concerned with the problem of efficient and accurate com- 
puter evaluation of (1.1) given w^(v) in tabular form, and trial 

values of 8. This problem is nontrivial, because, among other reasons, 
the operation must be repeatedly performed many times. 

In Section 2 we consider the class of quadrature formulae which 
would be exact if f(v,i0) were linear in the components of v. (The 
quadrature coefficients then depend on and are pre- calculated). 

A fairly sharp upper bound for the quadratures error is given, which 
depends (for given w. (v)} on the number and location of the quadrature 

points and the largest possible value of the mixed second partial 
derivatives of f(v'i8) with respect to the components of v. 
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If the instrument response windows are much longer in the 3 
direction than in either the v or a direction (which is apparently 
the case), it makes sense to consider a linear array of quadrature 
points in this direction only. 

For this case a feasible procedure is given whereby, for any 
fixed number n, the optimal set of quadrature points may be found. 

The optimal set of points is that set for which the error bound of the 
quadratures error is the smallest possible, assuming that a linear 
array is adequate. Upon finding the optimal set for each n under 
consideration, the corresponding error bound may be determined to see 
if it'is within the system accuracy constraints. Then the set with 
smallest n meeting the system constraints is selected. 

The adequacy of the linear array assumption may be checked 
by inspection of the final optimal set as determined above. A two 
dimensional procedure is described, for use if the one dimensional 
array is not adequate. 

In a simple analytic example the number of quadrature points 
required to meet a fixed error bound was found to be reduced by about 2,$ 
30% by use of the procedure given, as compared to equally spaced 
points. 

Although the procedure given may appear to require a large 
amount of effort, the more time consuming steps (involved with 
organization and numerical integration of the calibration data,), will 
have to be done for any accurate numerical integration procedure. 
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Section 3 briefly describes a least-squares procedure for 
estimating the Helium fraction, which probably has the advantage of 
being more stable than letting the Helium fraction be a free para- 
meter in the general non-linear least squares parameter estimation 
program. 

Section 4 gives an algorithm for linear estimates of the 
plasma parameters. A draft of this section was delivered informally 
in August, 1970. 

Section 5 is conclusions, and Section 6 contains three 
appendices. 
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SECTION 2. 


OPTIMIZATION OF THE NUMERICAL QUADRATURE 

2.1. General Considerations in the Optimization of the 
Numerical Quadrature, Introduction. 

Let f(v,6) be the particle velocity vector distribution 
with parameters 6, 9 = (8j ,0 2 ,. . . ,0^). Let w^(v) be the instrument 

response to a unit number density cold beam with velocity vector v 

fh 

at the k choice of collector, energy step and instrument position 
relative to the solar-oriented ecliptic coordinate system. In the 
discussion below we will suppress the subscript k. Thus w(v) = Wj,(v) 

for fixed k. w(v) is determined in the laboratory, for selected 
values of v, and the functional form of f(v,8), as well as bounds 
on the possible values of 9 are known. The problem considered here 
is to develop an efficient quadrature procedure for evaluating 
numerically 

1(0) = fff f(v,0)w(v)dv . (2.1.1) 

It is highly desirable to do this with as few quadrature points as 
possible, since, 1(e) must be evaluated (for each k) for many 
different values of 0, each time the plasma parameters are estimated 
iteratively by non-linear least squares. 

Let 

i(e) = l f(x v ,e) ///., v)dv (2.1.2) 

v=l R ' 
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where n is the number of quadrature points, RpR 2 ,...,R n are n 

disjoint regions in Euclidean 3 space to be determined, whose union 
is S, S being the region for which w(v) > 0, and x v , 
v = l,2,...,n are n vectors to be determined, x v e R^. Equation 


(2.1.2) is known as the quadrature formula. Iff w(v)dv is calculated 

R v 

and stored In advance, so only given linear combinations of f(x v ,e), 
for v = l,2,...,n are calculated during the iterative estimation of 
0. Then 


1(e) - 1(e) 



v=l 


R 


w(v)[f(v,9)-f(x v ,0)]dv . 


(2.1.3) 


V = 


By expanding f(v,6) is a Taylor series about x v , and letting 
(v-|,v 2 ,v 3 ), x v = (x^,x^,xp, we have 

3 


f(v,e)-f(x v ,e) = l (vj-xp 


0=1 »xj 


+ J_I I a?f(v.e) 


21 i j av^v. 


(v r xp(v r xp (2.1.4) 


v=v„ 


where v* is some point between v and x v , depending on v. 
Thus 
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i(e) - i(e) 



I 3f(x / e) - -^/w(v)(v r xpdv 
3=1 3Xj R v ' 3 3 


+ I i /// w(v) l ! 


v=l 2! R v i 5 l j=l 3v.8Vj 


v=v. 


<V*P<Y x 3 )dv . 


(2.1.5) 


v - . 


‘If x. is always chosen so that 

J 


in,. 


w(v)dv 


x v = v 

X j 


Jff w(v)dv 


, j = 1,2,3 , 


( 2 . 1 . 6 ) 


then the first term in (2.1.5) is identically 0. In Appendix A we 
show that, for a Maxwellian distribution 


max 

v, e,i,j 


3 2 f(v,e) 


9v. 3V 


j 


M 


(2.1.7) 


where 


M = 


(27r) 3/2 A- ] A 2 X, 


X . 
min 


( 2 . 1 . 8 ) 


kT. 


1 , where k is Boltzman's 


where N is the number density, , X^ = m 
constant, m is the particle mass, and T. is temperature along the 
direction of the eigenvector of the pressure tensor corresponding to 
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the eigenvalue X^ . is in units of velocity. If f(v,6) is the sum 

of two Maxwellian distributions as in (3.2.5) then (2.1.8) holds with 
N replaced by N p +2N a , Np and being the proton and a-particle 

number densities, respectively. Also, letting || x j x^+x|+Xg, 




= 3 tlv-x v |i 2 


< 3 max i!?-n|f , v e R. (2.1.9) 

max [|£-ri|j is the maximum dimension of the region R . Therefore, 
S»rieR v v 

if x^ is chosen as in (2.1.6), then 

J 

| I(0)~I(e) j < |m l max H?-n|| 2 HTv/(v)dv . (2.1.10) 

v=! e,neR v R v 


The main idea of the preceeding discussion (Eqns. (2.1 .3)-(2.1 .6) and 
(2.1 .9) -(2.1 .10)) has been adapted from a memo of Mr. John Day, of 
Informatics. 

Assuming that w(v) has units of current per unit number density. 


then 1(6) - 1(e) has units of current, M has units of (number 

n 

density x velocity ) , max [|£-nl| 2 has units of velocity , and 


/// 


€,neR. 


v 


-1 3 

w(v)dv has units of (current x number density x velocity ). 
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We now. proceed under the assumption that accuracy standards on 

the quadrature can be established based on 1) the requirements for 

accuracy in the non-linear least squares program 2) the accuracy of 

the calibration data w(v) and 3) other system errors (including 

digitization, timing uncertainties and the inability to accurately 

recover the signal current at the output of the collectors). Thus, we 

n 

wish to find n, and {R } so that n is as small as possible 

v v=l 

subject to the constraint 

I I(e)-i(e) | < y M l max j|S-n ip/// w(v)dv < C. (2.1.11) 
v=l C,neR v R v “ 

where is given (in units of current). 

Thus, we must find the smallest feasible and practical n, and 

n n 

corresponding sets {R } , {x v } for which 

v v=l v=l 

I max ||5-n|| l 2 /// w(v)dv < C 2 - 1 1 C, . (2.1.12) 

v=l S,neR v R v 

The remainder of this Section is divided into four subsections. In 
Section 2.2 we derive some theoretical lower Uounds for 

l max ||£-n|i 2 //Jw(v)dv . (2.1.13) 

v=l s,neR v R v 

The purpose of this is twofold. 1. The results can be used to tell us 
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whether it is worthwhile to attempt to improve a given choice 
n 

of 1 {R } , that is, how close is a given choice to the best 

v v=l 

obtainable. 2. The conditions under which the theoretical lower 

bound is attained give us some qualitative criteria for a good set 
n 

{R > for given n. In Section 2.3 we discuss qualitatively some 
v v=l 

of the calibration data available during the period of this study. 
Based on the criteria developed in Section 2.2, it appears that a 
satisfactory solution may well be obtained by exploiting the fact that 
the instrument response windows are much narrower in the radial , 

r 

(increasing ~) and ,azimuthal (a) directions than the elevation 


(B) direction. 

In Section 2.4, the results of Section 2.3 are used to formulate 

n 

the problem of choosing the ,{R } as a much simpler problem, namely 

Vl 


as a problem in determining boundaries in one dimension, instead of 

three. For this simpler problem we given an algorithm for obtaining 

n 

starting guesses for the boundaries of the regions {R } which 

v v=l 


minimize 


l max lU-n|| 2 Iff w(v)dv (2.1.13) 

v=l £ ,T}eR R 

V V 


for any given n. The calibration data is used to define w(v) here. 
The use of Newton’s method to iteratively improve the choice of 
boundaries to minimize (2.1.13), and the convergence properties of 
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this method, are described. For each candidate value of n, then, the 
best or near best set of boundaries may be found by the method described. 
Then the smallest n for which the inequality (2.1.12) holds may then 
be selected. Section 2.5 discusses a check for determining whether 
the one-dimensional solution is adequate. A two dimensional procedure 
is discussed in the (unlikely) case that it is not. In Section 2.6 
a simple example is discussed in which th? roughly 30% fewer quadrature 
points are required for the best choice of points as compared to > 
equally spaced points, for the same accuracy. 

2.2. Theoretical Lower Bounds for the Quadrature Error 
and Properties of a Good Choice of Regions 

By use of the Holder and Jensen inequalities we first obtain 
a theoretical lower bound on 

l max w(v)dv (2.1.13) 

v=l 5 ,neR v R v 

over all possible choices of regions {R^} whose union covers S. 

By examining the conditions under which this theoretical lower bound 
is attained, we obtain clues as to the properties of a good choice 
of regions. 

The volume V of a sphere in 3 dimensional Euclidean space 
is given by 

v - f-d 3 , (2.2.1) 
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where d is the diameter. Since a sphere is the geometric figure 
which maximizes the ratio of the volume to the longest dimension, we 
always have 


max [i£-n|j i2 >_C 



( 2 . 2 . 2 ) 


wh§re C = C 3 = (|r) 2 ^ 3 = 1.54 and equality is obtained if and 

only if R v is a sphere. In general , S, the region over which 
w(v) > 0, cannot be partitioned into spheres. If S is partitioned 
by, say plane surfaces only, then we will have C = = 2 2 = 4, 

obtained for R v a cube, since then V^ 3 = length of a side = 2 
x longest dimension. 

Therefore, provided only planes are used to partition S, 

then 



max Hs-n || 2 /j(jw(v)dv > C 4 l 

5*ne R v “ Vl 


1®.- 

r2/3 

Iff w(v)dv 

K 


L r v j 


(2.2.3) 


with equality holding if and only if {R^> are all cubes. 

The famous Holder inequality says, for any a, g positive 
numbers with a+B = 1, and a(v), b(v) positive functions. 


/// 


a“(v)b 8 (v)dv 


v 


'/// 


a(v)dv 


a 


Iff 

- R v 


b(v)dv 


(2.2.4) 
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with equality holding if and only if a(v) is a constant times 
b(v). Letting a = |- , $ = , a(v)sl, b(v) = w(v), we have 


/// w 3 ' 5 (v)dv 
R 



3/5 

w(v)dvl 


(2.2.5) 


and, taking the 5/3 power of both sides, we have 


Ilf w 3 ' 5 


(v)dv 


- R v 


5/3 


- R v 


IILf 3 riff 


w(v)dv 


R v 


( 2 . 2 . 6 ) 


and 


n r 




v 


Iff dvf / 3 r /// 


- R v 


w(v)dv 


1 C 4 l 


n r 


v=l 


Iff w 3 ' 5 


5/3 


(v)dv 


- R v 


(2.2.7) 


with equality holding if and only if w(v) is constant on each R v . 
By Jensen's Inequality (see Appendix B), we have 


n r 


'JIM 5 


(v)dv 




~ 5/3 i i ffi^r 

n V=1 R 


> c 4 n 


, 2/3 


7 /J> 


-r5/3 


(v)dv 


( 2 . 2 . 8 ) 


L S 


with equality holding if and only, if 


fff w^(v)dv = fff w 3//5 (v)dv 
R. R 


all u,v . ( 2 . 2 . 9 ) 
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In summary, for any partition of S into n disjoint regions 


n 

{R ,} whose union is S, a universal lower bound for (2.1.13) is 
v v=l 


given by 


n 

J, 


max 

?,neR. 


lls-1 


Ilf 


w(v)dv > 


V 




■Iff 


w 


3/5 


l5/3 


(v)dv 


. ( 2 . 2 . 10 ) 


The lower bound is attained if and only if equality obtains in 
(2sr2.2), (2.2.7) and (2.2.8)) which means (restricting the partitions 
to planes, to simplify the discussion). 

i) All the are cubes (that is, the ratio of 
volume to (diameter) 3 is maximized 

ii) w(v) is constant on each R^ 

iii) The volume of is inversely proportional 
to w 3/5 (x v ), x v e R v . 

This lower bound cannot be obtained for arbitrary w(v) since, 
for example i) and ii) require w(v) to be constant on cubes. However, 
from the above we may make the following semi -quantitative statements 
about properties of a good partition. 

i) The ratio of the longest to the shortest dimension of 
each should not be too large 
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ii) The regions should be chosen so that w(v) 
does not vary much for v e . 

iii) The volume of the regions should be small where 
w(v) is large, and conversely. 


2.3 Analysis of the Calibration Data and Simplification of 
The Problem 

The following discussion assumes, for concreteness, that the 
calibration coordinate system is a spherical coordinate system. In 
a spherical system a cold beam with rectangular coordinates (v-j,v 2 ,Vg) 

has velocity r = (v^+v|+v|) ly/r2 , and angles of incidence a and 3 
where 


v^ = r cos 3 sin a ( 

v 2 = r sin 3 (2.3.1) 

Vg - r cos 3 cos a 


See Figure 2.1. In what follows the volume element dv-j dv 2 dVg 

is replaced by the spherical volume element r 2 cos 3 drd3da. If the 
calibration coordinate system is not exactly spherical, the appropriate 
volume element should be used, but the remainder of the discussion 
below is unaffected. 

Let w,(v-j., 3 v 2 ,v 3 ) be the response of the instrument to a cold 
beam of protons of unit number density and coordinates v-j,v 2 ,Vg. 

w(v-|,v 2 ,v 3 ) is actually measured by 1) Fixing r and 3 and 



varying a. 2) Fixing r and d and varying 3. 3) Choosing 

a new r and repeating 1) and 2). Thus, it will be convenient 
to perform some calculations involving calibration data i f n the 
calibration coordinate system. 

Figures 2.2 and 2.3 reproduce some calibration curves and a 
plot which were available while this study was being performed. Com- 
plete labels and scale factors were not all provided and the missing 
information is omitted. Some labels have been added according to the 
coordinate system assumed in Figure 2.1. 

Figure 2.4 shows a hypothetical calibration curve as it might 
appear if the incidence angles a, B are kept fixed and the calibration 
bean velocity is varied. Based on discussions with ARC personnel, it 
is assumed that the resolution of the instrument in the (E/q) 
direction, i.e. in the direction of increasing beam velocity, is very 
narrow compared to the resolution in the direction of changing a 
or B. 

That is, the width of the curve in Figure 2.4 is very narrow 
compared to that of Figure 2.2a 3 when both are visualized in a common 
rectangular system. 

Figure 2.5 depicts schematically a and B calibration curves 
transformed to a rectangular system, for fixed r. Inspection of 
Figures 2.3 and 2.4 and their visualization in a common rectangular 
coordinate system reveals that the region S in a rectangular system for 
which w is not zero is shaped roughly like a possibly flattened hot dog. 
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A cross section of this hot dog shaped region, for fixed r, is the 
shaded region of Figure 2.5. The full hot dog shaped region S is 
depicted schematically in Figure 2.6. This is the region that must 
be partitioned. Due to the highly elongated shape of this region, it 
is consistent with the criteria for a good partition at the end of 
Section 1. to consider partitions perpendicular to the long axis of 
S. There is considerable benefit to being able to describe the 
partition simply in the calibration coordinate system. To this end, 
consider slicing (i.e. partitioning) the hot dog S with cones of 
constant 3. If this is done, then a Region R^ consists of all 

points whose polar coordinates (r, a, 3) satisfy 


(r, a, 3) e S' 


M e < <Vi 


(2.3.2) 


where the cones determined by 6 = 3 V and 3 = 3 V+ -| are two 
boundaries of R . Figure 2.7 shows a cross section of the region S 
in the r, 3 plane, for a = 0, and a cross section of R^ with 


boundaries 3 = 3 V and 3 = 3 V+ -|. 

Then, if we abuse notation by letting w(r,a,3) be the ' 
instrument response to a beam with spherical coordinates r, a, 3, 
then 



w(v-j,v 2 ,v 3 )dv^dv 2 dv 3 


6 v+l 2tt » 

= f cos 3d3 f da f w(r,a,3)r 2 dr 
3 a=0 r=0 
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(2.3.2) 



Of course the limits r = 0, r = » and a = 0, a = 2-n could be 
replaced by the extreme values of r and a inside S. 

Define 


2ir oo 

h(3) = cos g J da J w(r,a»8)r 2 dr . (2.3.4) 

a=0 r=0 


We may rev/rite (2.3.3) as> 



R 


v 


w(Vj ,v 2J v 3 )dv-jdV2dV2 


p v+l 


J h(B)d8 


(2.3.5) 


Vie next suppose that the longest dimension of R v will always be in the 

3 direction. If S is to be partitioned into, say 8 or 9 regions 
then this is reasonable since the hot dog S is around say 8 or 9 
times as long as it is thick, loosely speaking. In this case, the 
maximum dimension of R y is given by 

max ||5-n|[ = r*iS v+1 -8 v [ (2.3.6) 

£»neR v 


where r* is some intermediate value of r in S. See Figure 2.7. 
With the restrictions that a) S is to be partitioned by cones of 
constant 8 and b) the longest dimension of R v is in the 8 
direction, we may write 
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I max Hc-nlpJl/wtv, ,y,,Vo)dv-,dv 9 dVo' < 
v=1 e,veR v R v 1 L J ' 


C 4 (r**) 2 Jj^ SP V+1 



h(8)d8 


(2.3.7) 


where r** is some intermediate value of r inside S, and where 
3-j and 8 n+ -| are the extreme 8 values for points in S (see 

Figure 2.7). 8 is in radians in (2.3.7). If h(8) and n are 

given, it is possible to give an algorithm for choosing $2»^3* , ** , ^n 

to minimize 


n V+1 

L IW’J* / h(s)d6 • 

v 


(2.3.8) 


We do this in Section 2.4. To conclude this section we discuss the 
determination of h(8) from calibration data. 

Let 


i{i(a,8) ~ / w(r,a,8)r 2 dr . 


(2.3.9) 


r=0 


Then 


2tt 

h(8) = cos 8 J ^(a,’8)da 

a=0 


(2.3.10) 
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is the area under the curve r 2 w(r,a,g), as a function of r, 
where the curve w(r,a,$) as a function of r is shown* in Figure 2.4. 

If the calibration points are taken at a sufficiently fine grid, 

.. 6 (1) „(2) e (n B> a _ a (1) 0t (2) . a (n a> 

P P jP y CC Ot jCX j 9 


0) (2) (nj 

r = r ,r »r , then 


e (i) l a 


h(3^') ~ cos 


3=1 k=l 


(2.3.11) 


(3+1) Jj-lh a Jk+1) (k-1) 


r— ) ( 


r' '-r' 
2 


where w|r^,a^,3^) 


are measured values of w, and a v , 


(n+i) (0) (n+1) 

a , r , and r are suitably chosen. 


We may write 


h(3 u; ) = cos (T u l 

j=l 


(1) .(g^-c .(i) 


Ip (a vj/ ,3 v,/ ) (2.3.12) 


where 


- f „p),g(J), S U)) (,(«)*« 

k”*1 


*A r (k+l). r (k-1) 


°° ■ \ 

f w(r,a^»6^ j r 2 dr . (2.3.13) 
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If the calibration points in the r direction are sparse* 
as in Figure 2.4» then the approximation in (2.3.13) will not be good, 
and care must be taken to get sufficiently accurate values of 

Now, suppose we have determined h(g) for 

(^ ) (2) ( n g) (i+1) (i) 

3=3 ,3 ,...,3 p where 8 U ~ ; - B u; is' much smaller 

than the difference between two possible region boundaries S^ + -j - 8 y . 

Say B^ 1+ ^ - 8^ ~ 1°. We will then proceed below as though 
h(B) =0 for 3 < 3-j, and 3 > 3 n+ -j and is known continuously for 

3i £ 3 £ 8 n4 q. The next section gives an algorithm for minimizing 
(2.3.8). 


2.4 Choosing an Optimal Set of {82,33 

The problem now is to choose 8 2 
81 < 82 < < 3 n < 8 n+ -] to minimize 

f ^2’^3 ^ <*Wl 


t • * • » 8 n ) 


»8 3 »...>8 n 

f (82 ,83, • 



satisfying 
..,8 n ) given by 

h(8)d8 . (2.4.1) 


By the Holder inequality of (2.2.4) with a= \ , 8 = -j we have 

I ^v+l \|/ 3 v+l \~ p v+l 1 

I / ( / h (B)d8 >/ [h(8 )] T d 8 . (2.4.2) 

3 v 
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Hence, upon cubing both sides of (2.4.2) and summing, we have 

n ^+1 n / e v + 1 1 \ 3 

l (B U+1 -BJ / > I J [h(B>] T d6 . (2.4.3) 

v = 1 \ e / 

V V 

Equality holds j,n (2.4.3) if and only if h(g) is a constant for 

M » < Vi • 


By Jensen's inequality, we have 
n / Bv+1 1 C ~ n & v+l 1 

l J Wen 3 dB > n ij / [h(B)] 3 de 

=l\ / “ "v=l { 


& v+l 1 

,3 






u f3. 


v 


(2.4.4) 


with equality holding if and only if the {g^} are chosen so that 


Vl 


n+1 


j [MBjfdS = 1 j [h(S)fdB . (2.4.5) 

B„ B, 


The (unique) values 3|,g|,..,,B* for which (2.4.5) is 
satisfied may be found graphically by plotting 


H(8) = 


g 1 

/ (h(B)f df!/ 


B 


n+1 


I 


1 

[ti(S)] 3 dB 


(2.4.6) 


and graphically finding (3* which satisfies 
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H(p*) = , v = 2,3,.. . ,n . (2.4.7) 

This graphical procedure is illustrated in Figure 2.8. 

Now the set p|,p|, 8* found this way does not, in 

general, minimize f(8 2 >**'»8 n ) since equality does not, in general 

hold in (2.4.3). However, it is shown in Appendix C that 



< 


3 ! | fri- , 5 | 


4 


max h'(e) 


(2.4.8) 


where, according* ’to (2.4.4 )~- 

n 2 

is an absolute lower bound for • • >8 n )* For the derivative 

h‘ not too large, and !8* +1 -8*l‘ t small (remember p is in radians 

here), (2.4.8) shows that the chosen * * ’^n 5 ^ ounc * by so ^ vin 9 

(2.4.7) will be close to the best possible set. The next step is to 
calculate f(8|,...,8*) and determine if 


n+1 


/ [h(3)] 3 dg 


I H x C 4 (r**) 2 f(e^,...,e*) < C, . (2.4.9) 

(Recall Equations (2.1.11) and (2.3.7)). If not, n may be increased. 
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or a second order improvement in the choice of the (B v > attempted. 

Given the set (G|». . . ,g*) of starting values determined as above, 

various iterative techniques can be used to minimize 

numerically. Most likely an adaptation of the Marquardt algorithm in 
use in the Plasma Probe non-linear least-squares estimation program 
can be made to do the job. At this level, the improvement possible 
will, in general, be only second order. For completeness, we describe 
a relatively simple procedure for improving the set 
(g|,G|,...,B*). It is Newton's method in n-1 dimensions. It will be 

be convenient to use vector notation to discuss Newton's method. Let 


t 


( ^2 » * • * » 3^ ) a f(3 2* 


,G n ) given by (2.4.1) and let 


g(G) = 


(2.4.10) 
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We have 


Mi- 

as* 



i = 2,3, 


9 2 f 


3 i+l 


4 (e i+1 -3 i _ 1 ) h (3 i ) + 2 j 


3 i-l 


h(6)d3 


+ 


h ■ {e i J 


(Mi-i 



(2.4.11) 


(2.4.12) 


h(g i ) 

. . ,n 


(2.4.13) 
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a 2 f 

33. 3B 


i+1 


S H1 


3 2 f 


a6 1+i ap i 


= -2 


(3 i+1 -e i ) (h(3 i+ i )+ h (3i )) + f h(p)d3 


®1 


1 = 2,3,. .. ,n-l 


3 2 f 


38 - 36 . .. 
l i+J 


= 0 for j t -U0,1 
i = 2,3,. . . ,n-l 


(2.4.14) 


Necessary and sufficient conditions that f(£**) be a 
(local) minimum are that g($**) = 0 and A(8**) be strictly positive 
definite. (A matrix A is strictly positive definite if xAx' > 0 for 
all vectors x.) Here we then have 


f(t) = f(|**) +(6-£**)g(|**) + 1(I-6**)A(|**)(|-^**) + 

higher order terms (2.4.15) 


and, if g(8**) = 0, then f(6)-f(fl**) >0 in the neighborhood of 

6**. A sufficient condition that A(jp be strictly positive 
definite is 


a 2 ftl) afftSl > 
33f 30j 


j 


i,j = 2,3, ...,n . (2.4.16) 


This condition holds here if 


h'(0 i ) (3^ 


e i-l J2 " (3 i+1 


-e^ 2 i=2,3 n (2.4.17) 
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is positive or sufficiently small negative. 

If the variation in h(3) for 3.j_-[ < (3 < 3.j + -j is not large 

compared to h(p) for £ 3,£ 3 i+ -j 4 y , or, if 
(3 i+1 -3i) 2 z (e.-p. i ) 2 , then 


h'^) _(3 r P i _ 1 ) 2 - (P^-P^ 2 


(2.4.18) 


is small compared to the- remainder of the contribution to 
3 2 f 


, that is, to 


33- 


3 i+l 


4 (8 i+1 -3i_ 1 )h(3i) + 2 J h(8)d3 

Pi-1 


(2.4.19) 


and can be ignored. In the computations below it may be appropriate to 
do this. If not, h'(3^) should be estimated from the calibration data, 
possibly by first smoothing h(8). 

To derive Newton's method expand g(£) to the first order in 
a Taylor series about £*. Then 


g(£) = g(f*) + a(H*)(!-3*)' 


(2.4.20) 


^ Note that this is independent of the units in which 3 is 
measured. 
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(2.4.21) 


Then, to get t** for which g($**) = 0, we have 
g(!**) a o: g(l*) + A(fc) (?**-!*) ' 

and £** is obtained from 

1**4* = - A -1 (e*)g (I*) . (2.4.22) 

This procedure may be repeated with $** as a new starting guess. 

If A {t) is strictly positive definite for all £ with 

e 2 < 8 3 < ... < 3 n » then f(£) is a convex function with a unique 

minimum, and Newton's method is guaranteed to converge to the minimum. 
If, whenever h(3) is increasing at fr. , we have )>(3^ + i~3^ ) 

and conversely, then A is strictly positive definite. This will be 
the usual case for the starting guess |* = (3jj »$§»•• * 5 ^) > since, 

loosely speaking, the method for choosing selects 3* + -|-3* 

inversely proportional to the average of h^(3) for 3e[f3*,3* + -j]* 

Thus we may expect that this procedure will converge rapidly to values 
of which minimize (2.4.1). 

So far, we have assumed that the longest dimension of the 
regions found this way is in the 3 direction. If this turns out to be 
the case then if (2.4.9) is satisfied, then (2.1.11) will be satisfied. 
If S is much wider in the a direction than 3**i~3^*> where 
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(3|*,$| s ** •>£**) are the solutions found above then 2 dimensional 
partitions should be considered. These are discussed in Section 2.5. 

2.5 Partitions in Two Dimensions 

A two dimensional partition procedure assumes that there is 
only one quadrature point in the r direction but n 0 points in the 

P 

3 direction and in the a direction. Suppose a set 

• • • »$n*) boundaries for the quadrature regions has been 

found as in Section 2.4, which meets the accuracy specifications in 
the sense that (2.4.9) is satisfied. Then there is no point in con- 
sidering any partitioning in the a-direction, provided that 
^v+1 ' 9 reater ^ an angular v/idth of S in the ct 

direction. 

If the resulting partitions in the 3 direction are half or 
one third of the width of S in the a direction, then the simplest 
procedure is to divide all the regions R v in half by a contour of 

constant a down the middle. In the unlikely circumstance that the 
solution partitions in the 3 direction are very small compared to 
the width of S in the a direction, it is recommended that the 
problem be turned around to solve for partitions (ot|*,ct|*,. . . ,a**) 
in one dimension the same as was done for 3. 
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That is, substitute g(a) for h(g), where g(a) is given by 

TT °° 

g(a) = J cosBdp J w(r,a, 3 )r 2 dr , (2.5.1) 

3 =-m r ~0 

analogous to (2.3.4), and proceed as before. Then let the final 2 
dimensional partitions R , y= l,2,...,n', v-l,2,...,n have 

boundaries determined by (a, 3 ) e v , i 'f a** 5 a < a*|-j , 

3 ** £ 3 < 3 **^ . An (overly generous) error bound for the quadratures 

is given by the error bound determined individually in either the a, 
or 3 direction, if either direction always is the direction of the 
longest dimension of the regions. 

It can be shown that if ${a,3) defined by 

i|>(a,3) = i|>(ct,B) cos 3 (2.5.2) 

can be factored into the form 

^(a>3) = ^-j (a)?{> 2 (3) (2.5.3) 

then the procedure for solving the problem separately in each dimension 
will give the partition (in 2 dimensions) minimizing (2,1.13). 

(Recall, for example, that the function f(x,y) =1, a-j £ x £ a 2 ’ 

b-j < y < bg, and zero elsewhere, is the product of two functions 

f-j(x) f 2 (y), where f^(x) = 1 , a-j _< x £ a 2 , zero elsewhere, and 
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f 2 (y) = 1, b-j <_y £b 2 , and zero elsewhere. A bivariate normal curve 
suitably oriented also has this factorization property.) 


2.6. Relative Efficiency of the Optimal Set of Quadrature Points 


In this section we give a crude example to illustrate the 
comparison of equally spaced quadrature points in 6, that is, 

»v 

(e 2 ,e 3 ,...,e n ), where 


V- 1 _ yd 

& n + l- p n 


( 2 . 6 . 1 ) 


versus the 3* chosen according to 


3* satisfies H(3*) = ~- 


( 2 . 6 . 2 ) 


where H(3) is given by (2.4.6). 

For simplicity, and without loss of generality, let 6-j = 0, 


3 n+1 =1. As the example, let h(3) = 3 P , [Note that the results 

below will be the same for h(g) = KQ P and [0,1] replaced by any 
interval .] 

Then 


T(32>60> • ♦ * >3p) ^ f( 62**33* ’ 


••er> >-T 

n 


6 n+l 1 

/ Ce p 3 3 de 




_1_ 1 

n 2 (1+f) 3 


(2.6.3) 
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where • •>££*) are the optimal set. 

From the definition (2.4.1) of f^sgg,...^), 

1 

f (3 25 3 3 s — »3 n ) = —• / h(3)dg = — |p^pj 

" 0 


For p = 3, the solution to (2.6.2) is 


K - <"r> 2 


s v ~ 2,3,. . . ,n 


and 


f(8|»...,6*) 


n 

ii 


1 


1 


( SL) Z - (^-) 2 

K n J v n ' 


JF 

/ 6 3 dg 

J /¥ 


■ -L 5 C z - b-0*) 


4n 3 v=l 


(2v-l) 


(2.6.4) 


(2.6.5) 


(2.6.6) 
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~ ~ 1 

For p = 3, f(B„,...,S ) = — — . Letting 
1 n 4n 2 


e(n) = 


f ( $2 > • * • ) 


(2.6.7) 


we have 


e(n) = 1 l (v 1/2 -(v-0 ,/2 ) 2 (2v-l) 

n V=1 


( 2 . 6 . 8 ) 


9(n) is tabulated below for n = 4,5,. ..,10 


n 

4 

5 

6 

7 

8 
9 

10 


6{n) 

.631 

.605 

.587 

.575 

.565 

.558 

.552 


The absolute lower bound on f(6**,... ,6**) for p = 3 


1 


is — — x X , whereas f(g*,...,6*) = -1— 6(n) . 

4n 2 ^ L n 4n 2 

So for this case, f(B£,...,8*) is very close to the lower 

bound on the error, which is half that for equally spaced points. Since 
the error bound is inversely proportional to n 2 , then the general 
conclusion, for this example, is that the number of quadrature points 
may be reduced, over the equally spaced case, by a factor of at least 

' / l 0 ( n ) jz * 71 * 
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SECTION 3 


ESTIMATES OF THE HELIUM FRACTION 
3.1. Introduction 

In this section we briefly present two relatively simple 
methods for estimating the Helium fraction. The methods are based on 
the use of the same non-linear least squares parameter estimation 
program for which the quadrature methods of Section 2 were developed. 

We only briefly discuss the relative accuracy of these methods as 

compared to other methods, since it will depend crucially on how well 

\ 

the instrument resolves the proton and a-particle distributions. 
However, a comparative accuracy study should not be hard to perform, 
given a complete non-linear least squares program with data simulation. 


3.2. Assumptions and Definitions 

The particle velocity vector distribution f (v,6 ) for 

r r 

protons is assumed to be of the form 


f p (v ’ e p J 


-i(v-li„)r D, V (v-p„)' 

2' r p' p Xp p v ^p 


(2tt) 


3/2 


n 

i=l 11 


(3.2.1) 


where N is the proton number density, ii is the proton bulk 

r r 

velocity vector, r_ is an orthogonal matrix and D ln is a diagonal 

P Ap 
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matrix with diagonal elements 



= 1,2,3, where 




, i = 1,2,3 


(3.2.2) 


k is Boltzman’s constant, T. is the proton temperature in the 

i.L 

direction determined by the i tn column of r , and is the proton 
mass. 

If the instrument calibration is performed with a proton beam, 
then the particle velocity vector distribution of a-particles looks 

to the instrument like f (v,0 ) given by 

or a 


f a< v > 9 a> 


2N 


a 


1 




; ; ^/2 3 

(27t) n ^7" 
i=l la 


(3.2.3) 


where N , p , T , D. , X. , T. and m are defined for a-particles 
a K a’ a Xa’ la la a r 

analogous to the same quantities for protons. Assume now that 

def def def 

u = p = p, r = r = r and 4T- = T- . Then X- = X. = X . , 
K p % p a ip la ip ia l 5 


since 4m = m . Let g(v»p,Z) be defined by 
p a 


-1 


g{v,p,s) - 


-(v-p)z~ (v-p 1 ) 


( 2 w ) 3/2 | Z | 1/2 


(3.2.4) 


with Z = FD^r ' /. Then the sum of the two particle velocity vector 
distributions looks to the instrument like 
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V V ’V +f J V,6 cP = N p g( v »^> s ) + 2N a g(v^.E) • 


(3.2.5) 


iL 

Let w^(v) be the instrument response at the K" n combination 

of collector, energy step and relative instrument position, 
k = 1 ,2, . . . ,n.^ Then the k^ 1 data point, or instrument response 
I k , due to an input (3.2.5) is modeled by 

CO 

! k * I k ( V N «” I>i:) = N P 

—CO 

oo 

+ 2N a fJJ\(v)g(v,/2u,Z)dv, (3.2.6) 

-CO 

k = 1 ,2,. . . ,n . 


3.3. Description of the Method 

Let h^ = hy (u,S), i = 1,2, k = l,2,...,n, be defined by 

h kl (vi,S) = iUw k (v)g(v,y,2)dv (3.3.1) 

k = l,2,S..,nw 

h k2 (y,s) = 2 JJJ w k (v)g(v,v^y,E)dv 


The system of equations (3.2.6) may be rewritten 


3.3 

^ This n is unrelated to n of Section 2, we are running out of symbols. 





(3.3.2) 


3.3.1. Equations for Method 1 

Let H « H(p,z) be the n x 2 matrix with ki^ entry h^. . 
If y,E are assumed, and the data 1^ ». . . , I are given, then, using 

A A 

(3.3.2), the least squares estimates N and N for N and N 

papa 

are given by 


M 

V 

A 


1 'A 
• 

* 

= (H , H)“ 1 H* 

N 


i , 

7 


7 


The Helium fraction 6 is estimated by 


H = H(y,2) . 


N /(N +N ). 
a ' p a 


(3.3.3) 


3.3.2. Equation for Method 2 « 

If N, y,E and I^,...,I n are given, the least squares 

A 

estimate 6 for <5 is, from (3.3.2) 
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1 


(3.3.4) 


6 


i | 1 ^ h i2~ h i1 )" 


N J 1 (I i“ h il )(h i2~ h il 


3.3.3. The Procedures 

Procedure for Method 1: Obtain starting values for p and £ 

by the use of the non-linear least squares program with N q assumed 
0. [Delete all obvious 100% a-particle data from the data set.] 

A A “* , 

Using these p and £, determine N p and N a via (3.3.3). introduce 
the model 

f p (v,e p )+f a (v,0 a ) = N p g(v,p,£) + N a *2g(v,/2p,E) (3.3.5) 

A A 

into the non-linear least squares program with N and N fixed, 

p ot 

and p, 2 as parameters to be fitted. Using the updated estimates of 

A A 

p and 2, obtain new N and N from (3.3.3) ; . 

p 0» 

Procedure for Method 2: Obtain starting guesses for N, p, 

and E as in Method 1. (Starting guess for N = starting guess for 
N p ). Using these values for N, p and £ estimate 6 by the use of 

A A A 

(3.3.4). Iterate, if necessary by replacing N and N by N(l-S) 

p ex 

A A 

and N<5 in (3.2.5) with N a free parameter and 6 fixed. 
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3.4. Conclusion 


Either method is simple, assuming the existence of auxiliary 
programs which must be developed in any case. The general model (3.2.5) 
with all parameters to be determined may also, of course, be used in 
the non-linear least squares program. In the case where the two 
distributions are poorly resolved by the instrument, then there is some 
reason to believe that either of the above methods, particularly 
Method 2 if 5 is small, v/ould tend to be more stable than fitting 
the full general model simultaneously. 

For the case of the two distributions resolved to the extent 
that two distinct peaks are identifiable, these methods should be 
compared with methods based on the relative peak heights, etc. For 
the case of the distributions completely resolved, may be 

A 

estimated in a least squares sense by N given by 

Cl 


A 





(3.4.1) 


where the summation is taken only over k for which 1^ contains 

no contribution from protons, and y and 2 are determined 
independently from the proton data. 


3.6 



SECTION 4 


LINEAR ESTIMATES OF PLASMA PARAMETERS 


4.1. Introduction 

Linear estimates of the plasma parameters (for data from protons 
only) can be made, because the parameters are simple functions of the 
moments, which can be viewed as linear functionals. The linear 
estimates can be expected not to be as accurate as the non-linear 
estimates when the Maxwellian distribution holds, since the latter 
uses a known functional form for the particle velocity vector distri- 
bution. However, the linear estimates will be cheap to compute, since 
the coefficients can be calculated ahead of time. Also the linear 
estimates of the moments are reasonable estimates of the moments 
irrespective of the true nature of the particle velocity vector dis- 
tribution. Their accuracy (not counting the errors induced by poor 
calibration data) may be examined ahead of time by simulation. 

Let f(v,eKv' be the particle velocity vector distribution, 
where 0 stands for the parameter set 0 = {N, y,Z}, as defined in 
previous sections and v = (v-|»v 2 >v 3 ). (We are assuming data from 

protons only). 


If 


f(v,e) 


N - i(v-y)E _1 (v-u ).* 

( 2*) 3 ' 2 | I | 6 


(4.1.1) 
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then 


p. 


cr. • 
i J 


Ilf f(v,6)dv 

siF v i f(v ’ e)dv ■ 1 * i > 2 ' 3 = 

i/JJ v i v j f(v ' 9 > dv = 


= nuinber density 

coordinates of bulk (4.1.2) 
velocity vector 

entries of the 
pressure tensor, 

2 = {a 13 > . 


As before, let w^Cv), k = l,2,...,n, be the current at the 

output of the target at the k observation step due to a cold proton 
beam of unit number density and velocity vector v. The observed data 
points 1^, k = l,2,...,n, (that is, current at the output of the 

target) are given by 


\ s Iff w k (v)f 0 (v)dv , k = l,2,...,n. (4.1.3) 

It will be convenient to use the inner produce notation (f,g), 

(f,g) =/// f(v)g(v)dv . (4.1.4) 

Let 

g 0 (v) - i 
g.j(v) = v i 

9 ij (w) * Vj 


, i = 1,2,3 (4.1.5) 

» i,j = 1,2,3 . 
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Then, writing f 0 for the function of v defined by f(v,0). 

N - (S o .f 0 ) 

Nm, = (g t ,f 9 ) i= 1,2.3 (4.1.6) 

We will find linear least squares estimates of the 10 quantities 

Wj , i = 1 , 2,3 5 N (cr.j ) s i <. j s 1 » 2 , 3 . 

Estimates of the parameters N, p. and a., can then be computed 
from these quantities. 

4.2. The Method in One Dimension 

To help in understanding the method as well as its pitfalls, 
we describe it in some detail for functions defined on the (one- 
dimensional) real line, as opposed to 3-dimensional velocity vector 
space. 

In one dimension we wish to estimate 


n = /g 0 (v)f(v,e)dv 

g 0 (v) s i 


pN =/ g-|(v)f(v,0)dv 

g 1 (v) = v 

(4.2.1) 

N(c 2 +p 2 ) = / g 2 (v)f(v,e)dv 

g 2 (v) = v 2 . 
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See Figure 4.1. In order to carry out the method we must have a finite 
region of integration outside of which f(v,0) is known to be always 
negligible. Call this region R. We will let 

(f»g) = / f(v)g(v)dv . (4.2.2) 

R 


Let ^ n be the (assumed n -dimensional) linear space of all 

n 

finite linear combinations of the functions {w b (v)> . Let us 

K k=l 

A 

comceptually approximate g, (v) by an element g. (v) in the space 

(i = 0,1,2). 


g i (v) = J, c ik w k (v) 


(4.2.3) 


where the {c.^} are to be found. Choose the {c.^} to minimize 


[jg-i (v) - g i (v) jp 


/ ( g i (v) " J-| c ik w k (v >) dv 

R (4.2.4) 

/ 3,00- 2 ^ c 1k a. k * g c 1k c u b kJ , 


where 


a ik 



/ 9 1 (v)w k (v)dv i = 0,1 ,2 

R (4-2.5) 

/ w k^ w £^ dv • 

R 
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The coefficients Cc^-[,c^ 2 >**-»c. n ) which minimize (4.2.4) are 
easily shown to be given by 





( a il * a f 2* — * a in) B 


-1 


(4.2.6) 


where B is the n x n symmetric matrix with k£ entry b^. 

B will be diagonal if the instrument "windows" do not overlap (in 

A ^ 

velocity vector space). Note that g^(v) will approximate g^ (v) 

very nicely in the region of interest if the {w^v)} are as in 

Figure 4.1 but not so well if they are as in Figure 4.2. See 
Figures 4.3 and 4.4. 

Once the coefficients (c 1 *i»c. 2 s...»c- ) , i = 0,1,2, are 
chosen, we approximate the desired moment 


/g i (v)f(v,0)dv (4.2.7) 

R 


by 


/ g,(v)f(v,e)dv = c jk /w k !v)f(v,9)dv = c ik I k 


i = 0,1,2 . 


, (4.2.8) 


An error bound is given by the formula 



/(g 1 (v)-g i (v))f(v,e)dv 


R 
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(Cauchy-Schwartz Inequality) 

o r - 1 1 /2 _ - 1/2 1/2 

j r (g^Cv)-g 1 -(v)) z dv. J f 2 (v,e)dv = |[g 1 -(v)-g 1 -(v)|| ||f(v,e}|| 

R R J 

(4.2.9) 

Note that the function 

f(v,e) = l f k w.(x) (4.2.10) 

k=l K K 

where 

(f r f 2 f„) = (i,(e) i n (e))B _1 (4-2.il) 

is an approximation to f(v,e) in «. . ^ n and satisfies 

(g i (v)~g i (v) , f(v,6)) * 0 . (4.2.12) 

Hence we may replace the right hand side of (4.2.9) by 

A 

/(g i (v)-g i (v))(f(v,e)-f(v,e))dv 

R 

5 l|g 1 -(v)-g i (v)|| 1/2 (|f(v,e)-f(v,e)|| 1/2 . 

Of course ||f(v,e)-f(v,0) || 1y/2 is not generally known in advance. 

Thus the error bound for the moment estimates depend on how well g^v) 
and f (v,9) may be approximated by linear combinations of the window 

(i 

functions {w,.(v)} 

k i._i 
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4.3. The Estimates in Velocity Vector Space 

The formulas are unchanged for three dimensional space. It 
is necessary to establish numerically from the calibration data the 
quantities 



b k£ = /yy w k (v)w Jl (v)d v . (4.3.2) 

R 

These are independent of observed data, of course, and are only 
calculated once. The coefficients (c Q -| ’ c 02’ * * * ,c on^ »* * * 

( c 21 ’ c 22 , ^ ,, * ,,c 2n’'^) usec * ™ ^titrating the moments, analogous 

to (4.2.7) and (4.2.8) are given by (4.2.6). 


4.4. Conclusions 

Clearly the accuracy of the method depends on the particular 
instrument design as well as the unknown particle velocity vector 
distribution. Errors in the calibration data and in the quantities 
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(4,3.1) and (4.3.2) will also affect the results. 


It is recommended that a pilot check of the accuracy of the 
method be made as follows. A set of test cases f(v,0) are chosen, 
(most reasonably Maxwellian). To study errors due to the method (as 
opposed to errors due to faulty calibration data, the w }< .(v) may be 

defined by smoothed calibration data. Using an accurate quadrature 
procedure, the quantities in (4.3.1) and (4.3.2) are calculated. The 
"data" 1^ should be simulated via the same calibration data and 

quadrature procedure 

h = //J w k (v) f(v,e)dv . (4.4.1) 

The moments and parameters are then estimated via (4.1.5), 
(4.1.6) and the appropriate 3 dimensional analogue of (4.2.7) and 
(4.2.8) obtained by erasing the in g^ and g. . 
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SECTION 5 


CONCLUSIONS 


1. Optimization of Numerical Quadrature. Probably the most 
difficult and time consuming part of the development of any accurate 
quadrature procedure using a minimal number of quadrature points is the 
determination of the volumes JJJ w(v)dv or JJJ w(v)dv, of the 


R v S 

calibration curves. Once this work is organized on the computer then 
the development of optimal quadrature formulae as described in 
Section 2 is not much more if any additional effort. If a quadrature- 


_ routine is to be u§|d for standard production work over a long period . 
of time, it is probably worth the investment to optimize this routine. 


The discussion of this section has illustrated an appropriate 
and relatively simple way to use the calibration data in the 
quadrature routine. 

2. An Alternative Procedure is described for estimating 

the Helium fraction. It will require some simulation to compare this 
procedure against other candidates. It is probably more stable than 
estimating the Helium fraction simultaneously with the other parameters 
in the non-linear least squares program. 

3. A method for obtaining estimates of the plasma parameters 
without iteration has been given. While this method is no doubt 
cheaper than iterative non-linear least squares methods, its accuracy 
has not been established here. It may be used, however, to provide 
starting guesses for an iterative method. Furthermore, it provides 
reasonable estimates of the moments of the distribution irrespective 
of its functional form, and thus may provide some descriptive numbers 
when the distribution is not Maxwellian. 
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SECTION 6 


APPENDICES 
APPENDIX A 

Upper Bounds on the Mixed Partial Derivatives 
of The Maxwellian Density 


Let 


f(v,e) = 




- ^(v-ji) _1 (v-y}) < 


A.l 


where 8 = (N,y,2), v = (v-|,v 2 ,v 3 ), a velocity vector. We show that 


where 


max 

^1 ,V 2 sV 3 , ®’i 


M = 


3 2 f 


3V • 3V . 
* J 


< M 


A. 2 


(2tt) 3/2 n AT 
i=l 1 


min 


A. 3 


where is the smallest among (X-pXgjXg), the eigenvalues of the 

pressure tensor 2. 


Represent 2 as rD^r 1 where r is orthogonal and D^ 
th kT i 

is diagonal with ii xn entry X. = . Let 


gU) = e 


l 3 x i 


, x = (X-j.Xg.Xg) . 4.4 
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Then 


f( v ,9) = L g(x) A. 5 

M 3 ' 2 n AT 

1=1 1 

with 

x = (v-y)T A. 6 

(IC = |E| 1/2 ). 
i=l 1 

Since the Jacobian of the above transformation from v to x is 1 
then 


max 


3 2 f 

= max 

3 2 q 

3V.3V. 

■ J | 

3x- 3x . 

* J 



(2t) 3/2 n 7 \. 

1=1 1 


A. 7 


The maximum curvature of g occurs at x i = x 2 = x 3 > in 


the direction of the x^ axis where is the smallest among 
^ 1 s ^2 s ^ s and so 


since 


max 

i » J »x. jXj 


. 9 2 g 

= sfa 

3x.3x. 

' J 

9x k 


xi=x 2 = x 3=0 




A. 8 


A. 9 


It is obvious that if f(v,6) is the sum of 2 Maxwellian distributions, 

as in (3.2.5), then A. 3 holds with N replaced by N +2N . 

pa 
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APPENDIX B 


Jensen's Inequality 


Let g(x) be any strictly convex function of the real variable 
x, and let Xj,x 2 ,...>x n be any given n real numbers in the 

domain of g, with 

n 


Then, 


5 " "ii • 


l m SU V ) > ng(x) 


B.l 


v=l 


B.2 


with equality holding if and only if x v = x, v = l,2,...,n . 

Proof : 

The straight line y(x) given by 

y(x) = g‘ (x)(x-x)+g(x) B.3 

is tangent to the strictly convex function g(x) at the point x = x 
and is below g(x) otherwise. Thus 

n n 

I g(x ) > g‘(x) l (x -x)+ ng(x) = ng(x) 
v=l v ” v=l v 

with equality holding if and only if x v = x, v = l,2,...,n. 
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To apply Jensen's inequality to (2.2.8), set g(x) 


ffjvt 3 / 5 (v)dv, and to apply to (2.4.4) set g(x) 

3 V+1 

/ [h(fi)] 1/3 dfS . 



APPENDIX C 


Proof of the Inequality in (2.4.8) 


We show that 


6 ^i 


I (^,-B*) 2 f He ) ds - - 

v=l n 


8* 

V 


Vl 

/ [h(e)] 1/3 ds 


£ I |8J +1 -6*I 4 mx |h'(e)| 


C.l 


where the g* satisfy 

ej-n 

/ [h(8)] 1/3 de = 1 / [h(e)] 1/3 de . c .2 

8* g, 

H v M 1 

Proof: 

If the {g*} satisfy the above condition, then, according to 
Jensen's Inequality 

C.3 

where we are writing g^ = g^ , g* +1 = 8 n+ ] • 

Thus, it is only necessary to show that 
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(6^,-e*) 2 / h(e)de - 


ft* 




/ [h(6)]^ 3 dB 


3* 


C.4 


max ! h ' ( 3) i (B^i-3*) 4 . 

6 ^ v ,fW 


By the mean value theorem, there exists a 3 V between 


B* + i such that 


ejti 


hl/3 (e v ) 58 / CK3)] 1/3 dg 


3* 

v 


and so 




B $+i 


/ Dite )] 1/3 de 


3* 

v 


Now 


h(3) = h(3 v ) + (3-3 v ) h'U) 


3* and 
v 


C.5 


C.6 


C.7 


where 5 e(3 v ,3 v+1 ) » and so 

|h(e) - h(5 )| < |6* *B | max h' (0) . C.8 

v - V+1 Vo ^,) 

Substituting (C.6), (C.7) and (C.8) into (C.4) gives the 


result. 
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